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Abstract. This work is part of an ongoing research programme to study possible 
Primordial Black Hole (PBH) formation during the radiation dominated era of 
the early universe. Working within spherical symmetry, we specify an initial 
configuration in terms of a curvature profile, which represents initial conditions for 
the large amplitude metric perturbations, away from the homogeneous Friedmann 
Robertson Walker model, which are required for PBH formation. Using an 
asymptotic quasi-homogeneous solution, we relate the curvature profile with the 
density and velocity fields, which at an early enough time, when the length 
scale of the configuration is much larger than the cosmological horizon, can be 
treated as small perturbations of the background values. We present general 
analytic solutions for the density and velocity profiles. These solutions enable us 
to consider in a self-consistent way the formation of PBHs in a wide variety of 
cosmological situations with the cosmological fluid being treated as an arbitrary 
mixture of different components with different equations of state. We show that 
the analytical solutions for the density and velocity profiles as functions of the 
initial time are pure growing modes. We then use two different parametrisation 
for the curvature profile and follow numerically the evolution of a range of initial 
configuration. 
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1. Introduction 

The possible existence of primordial black holes (PBHs) was first proposed in 1966 by 
Zeldovich & Novikov pQ and, at the beginning of 1970s, by Hawking 0. It was then 
widely discussed in the following years, see for example the recent review by Carr for 
a full list of references In 1974 Hawking made his famous discovery of black hole 
evaporation ^ , that is cosmologically relevant if the mass of the black holes concerned 
is less than 10 15 grams. For this reason the problem of PBH formation started to be 
attractive and it was widely investigated in the following 30 years. 

The PBH formation process was first investigated by Carr (1975) [5] using a 
simplified model for an overdense collapsing region, described as a closed Friedmann- 
Robertson- Walker (FRW) universe, surrounded by a spatially flat FRW expanding 
background. In the radiation dominated epoch of the Universe this leads to a threshold 
value for the perturbation amplitude S c , where the amplitude 5 is defined as the mass 
excess in the overdense region, and the black holes formed have masses of the order 
of the horizon mass at the time of formation. In this way Carr obtained a first rough 
estimate of 6 C ~ 1 /3 (evaluated at the time of horizon crossing, when the overdensity 
region enters into the cosmological horizon) by comparing the Jeans length with the 
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cosmological horizon scale at the time of black hole formation. Subsequently a self 
consistent hydrodynamical analysis of PBH formation was carried out by Nadezhin, 
Novikov and Polnarev in 1978 [5] and in 1980 [7] using, for the first time in this context, 
a hydrodynamical computer code written in the Misner-Sharp slicing, characterised 
by a diagonal metric with a cosmic time coordinate that reduces to the FRW metric 
in the absence of perturbations. Previously the same slicing was used by Podurets 
[H] and May & White to study stellar core collapse. An alternative analysis 
was developed in 1979 by Bicknell & Henriksen JH] using a different method based 
on integration along hydrodynamical characteristics. These papers showed that the 
threshold amplitude of the perturbation is dependent on the particular shape of the 
initial conditions. It was also found in [H] that pressure gradients considerably reduce 
the PBH mass formed at the end of the hydrodynamical process. 

In the following twenty years attention was concentrated on different aspects of 
PBH formation. For example, calculating the amount of Hawking radiation emitted 
by sufficiently small PBHs, with a mass less than 10 15 g, one can obtain important 
constraints on parameters that characterise the different epochs and processes of the 
Universe (see, for example, ^JE] and references there). The formation of PBHs was 
considered also within different scenarios, for example during phase transitions |13| . 
with a soft equation of state by collapse of cosmic loops ^H], JB! or from bubble 
collisions [T7J. In general, the study of PBHs provides a unique probe for different 
areas of physics: the early Universe, quantum gravity, gravitational collapse and high 
energy physics. This is explained with full lists of references in the various reviews of 
PBHs, sec for example 0. 

More recently, in 1999, Niemeyer and Jedamzik |18l IB?] made new numerical 
calculations pointing out the relevance of scaling laws for PBH formation. They 
showed that the black hole mass A/bh follows a power law (S — S c ) 7 if 5 is close 
enough to 5 C , the same behaviour as seen in critical collapse by Choptuik [201 and 
other authors (see review [23)- Niemeyer and Jedamzik found that S c ~ 0.7 for the 
three types of perturbation profile that they studied. In the same year, Shibata and 
Sasaki [22] presented an alternative formalism for studying PBH formation focusing 
on metric perturbations, rather than density perturbations, as had also been done 
previously [S] |23I - They pointed out that the initial conditions used in [Fij\ were 
specified initially within a nonlinear regime of perturbations of the energy density and 
velocity field, and were therefore "inevitably contaminated by an unrealistic decaying 
mode" component that would diverge for t — ► 0. In a recent paper [21] this analysis has 
been converted in terms of the perturbation amplitude 5, showing that the Shibata and 
Sasaki results, corresponding to a wide choice of perturbation shapes, are consistent 
with 5 C in the range 0.3 < S c < 0.5. 

The disagreement between this and the value S c ~ 0.7 has been explained by 
Musco et al [21], where simulations similar to those used in JH] were carried out but 
specifying the initial conditions within the linear regime, giving only growing mode 
solutions at the horizon crossing time. The simulations in [21] give values of S c in the 
range 0.43 — 0.47 (instead of 0.67 — 0.71) for the same types of perturbation profiles 
used in [T5| . 

The present work is a new analysis of PBH formation using the numerical 
technique developed in [23], implemented together with a quasi homogeneous solution 
( |30|. [5]) that, assuming spherical symmetry, can be characterised by a single function 
of the radial coordinate. In this paper this function, called K(r), is chosen as a 
curvature profile that allows a self consistent determination to be made of the whole 
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set of initial conditions. This is a further development of the original approach of 
Nadezhin, Novikov and Polnarev 0. 

According to the asymptotic quasi homogeneous solution ,30 (as t — > 0), an 
arbitrary curvature profile K (r) corresponding to a large amplitude perturbation of 
the metric does not depend on time, while perturbations of energy density and velocity 
vanish asymptotically as t — > 0. Therefore these can be treated as small perturbations. 
Solving the equations for the small velocity and density perturbations analytically 
we impose in a self consistent way all of the initial conditions corresponding to a 
moment in time when the quasi homogeneous solution of a certain order is valid. The 
curvature profile K(r) appears as the source in the right hand side of the relevant 
equations and we can say that the density and velocity perturbations which we use 
as initial conditions are generated by the curvature K(r). Then we use the computer 
code to follow the subsequent non linear evolution of the initial configuration (by the 
configuration we mean the region of strong metric perturbation) . To avoid confusion 
we should emphasise that the small perturbations predicted by any cosmological model 
are relevant in this context only when one calculates the probability of finding a 
configuration with a high amplitude perturbation of the metric. In fact, the small 
perturbations of density and velocity which we discuss here are determined not by 
arbitrary cosmological initial conditions, but entirely from the curvature profile within 
the initial configuration. 

The paper is organised as follows. 

In section 2 we give a very brief description of the Misner-Sharp equations. We 
then use these equations to specify all of the initial conditions, which we then use in 
our numerical computations, in terms of a curvature profile. When the length scale 
of the configuration is much larger than the cosmological horizon, the Misner-Sharp 
equations can be reduced to a system of linear differential equations which we solve 
analytically. 

In section 3 we obtain solutions for a rather general multi-fluid case, where the 
equation of state corresponds to an arbitrary mixture of perfect fluids. We demonstrate 
that, in the simple case of single perfect fluid, corresponding perturbations of energy 
density and velocity, obtained within the framework of the quasi homogeneous 
solution, evolve with time like a pure growing mode in standard cosmological 
perturbation theory, while their space dependence is entirely determined by the 
curvature profile K{r). 

In section 4 we discuss the physical properties of the curvature profile K(r), 
which is directly connected to the time-independent comoving curvature perturbation 
1Z frequently used in literature. 

In section 5 we introduce two different parametrisation of K(r). 

In section 6, we present numerical tests to demonstrate self-consistency of 
the initial conditions. We then show that the initial conditions described by the 
quasi homogeneous solution have been imposed consistently in the code and give 
numerical examples of primordial black hole formation in the case of a radiation 
dominated universe. We show how the threshold for black hole formation is linked 
with the curvature profile and discuss the results obtained with the two different 
parametrisation used in the computations. 

Summary and conclusions are given in section 7. 
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2. Mathematical formulation of the problem 

2.1. The Misner- Sharp equations 

Assuming spherical symmetry, it is convenient to divide the collapsing matter into 
a system of concentric spherical shells and to label each shell with a Lagrangian co- 
moving radial coordinate which we denote as r. Then the metric can be written in 
the form used by Misner & Sharp 120] 

ds 2 = -a 2 dt 2 + b 2 dr 2 + R 2 {dO 2 + sin 2 0dip 2 ) , (1) 

where R is a circumference coordinate, a and b are functions of r and of the time 
coordinate t , which reduces to the familiar FRW time coordinate (referred to in 
literature as "cosmic time") in the absence of perturbations. The FRW metric used 
to describe homogeneous and isotropic cosmological models is a particular case of Q : 

dr 2 



ds 2 = ~dt 2 + s 2 {t) 



r 2 (dO 2 +sin 2 



where s(t) is the scale factor and K is the curvature parameter that is equal to 0, +1 
and —1 for flat, closed and open universes. 

For a classical fluid, composed of particles with nonzero rest-mass, it is convenient 
to use the rest-mass p, contained interior to the surface of a shell (or, equivalently, 
the baryon number) as its co-moving coordinate r. For fluids not possessing these 
conserved quantities, one can still define a "relative compression factor" p |271 I28) 
(equivalent to the rest-mass density on a non relativistic fluid) and one then has 

dfi = 4irpR 2 bdr . (3) 

Identifying /i and r one has 

b = iiy (4) 

Following the notation of |26) . we write the equations in terms of the operators 



m 

b \d/j> 



and define 



U = D t R, (7) 

T = D r R, (8) 

where U is the radial component of the four- velocity in the associated Eulerian frame, 
R is the circumference coordinate, and T is a generalization of the Lorentz factor. 

We assume that the matter consists of different components and that each 
component can be described as a perfect fluid with the equation of state 

P = ie, (9) 

where p is the pressure and e is the energy density. For any equation of state of the 
form p = p(e), the system of Einstein equations can be written as: 



D t U 



M 

-D r p + — + AirRp 



(e+p) rl ~ ' i? 2 



(10) 
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Dtp 



P 



D r (R 2 U), 



(11) 



TR 2 



D t e 



e+p 



Dtp, 



(12) 



P 



D t M = inpUR 2 , 



(13) 



a 



(14) 



D r a 



e + p 



DrP, 



D r M = AnTeR 2 , 



(15) 



where M is a measure of the mass-energy contained inside radius /i and T can be 
calculated either from JSJl or from the constraint equation 



When the cosmic time approach is used for calculations of gravitational collapse 
leading to black hole formation it has the well-known drawback that singularities 
appear after a finite time in the calculations, before an event horizon has formed. 
Once a singularity has formed, the calculations cannot be continued further and so 
parts of the evolution which could potentially be seen by an outside observer cannot 
be followed in this way. In particular, it is not possible to follow all of the process of 
the formation of the event horizon. This drawback is the reason for an observer-time 
formulation being used in the present work rather than continuing to use the cosmic 
time formulation. Hernandez & Misner 29 introduced the concept of "observer time" , 
using as the time coordinate the time at which an outgoing radial light ray emanating 
from an event reaches a distant observer. In the original formulation, this observer 
was placed at future null infinity but for calculations in an expanding cosmological 
background we use an FRW fundamental observer sufficiently far from the perturbed 
region so as to be unaffected by the perturbation. A complete description of the 
Hernandez-Misner equations was given in |25| . 

2.2. Perturbations in the quasi homogeneous solution 

To describe PBH formation, we follow the quasi homogeneous approach ([H], [7]) and 
solve the Misner-Sharp system of equations, giving a non perturbative description 
of large amplitude metric perturbations away from the homogeneous FRW model. 
Using this approach, as mentioned in the introduction, we can avoid the unphysical 
arbitrariness in the choice of initial conditions considered in previous works |18U22ll25| . 

The characteristic feature of the asymptotic quasi homogeneous solution is that 
all mass elements expand asymptotically, when t — > 0, according to the FRW model, 
with energy-density e = 1/6(1 + j) 2 irt for an equation of state given by (0, while a 
spatial hypersurface with t = const can have arbitrary curvature if the perturbation 
has a length-scale sufficiently larger than the cosmological horizon. As a consequence 
of this the hydrodynamical quantities can be considered as small perturbations with 
respect to the background solution, while the curvature perturbation is arbitrary large. 
This curvature perturbation is time independent, because pressure gradients in this 
regime are negligible ([23], M see also, for example, [2]). 

In this work, we consider spherical symmetry and specify conditions on an initial 
space-like hypersurface using a time independent function K = K{r) that represents 
the curvature profile. As we will see in the next section, this K(r) is directly related 




(16) 
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to the comoving curvature perturbation 1Z. We can also see that introducing a profile 
K{r) into the Friedmann Robertson Walker metric one obtains an asymptotic 
solution of the Einstein equation in the limit t — > 0. We therefore use K (r) to represent 
initial curvature perturbations with scales much larger than the horizon, and solve 
the system of perturbed differential equations to express the set of hydrodynamical 
variables in terms of K{r). 

The system of equations (0 - (|lf>|) can be re-written as: 

R = aU (17) 

- = -T 2 -- (W) 

a 1 + 7 e v ; 

M = -4:TTjeR 2 R (20) 
M' = AireR 2 R' (21) 

where the dot and dash denote differentiation with respect to t and r respectively, 
we use the equation of state Q to express the pressure as a function of the energy 
density, equation |JSJ gives an expression for T, and p can be calculated from equation 

CD- 

The background solution is a spatially flat FRW universe described by K = 0. The 
corresponding value of the energy density (indicated with the suffix "b" ) is calculated 
from the Friedmann equation 

The quantity e b is related to the scale factor s(t) by equation (|2U)l used for the 
unperturbed case 

^ = -3(1 + 7)", (24) 
e b s 

while the background values of the other quantities are obtained from the FRW metric 
JSJ and the set of equations l(T7jl - (|2~2"1) : 

a h = I (25) 

b h = s{t) (26) 

Rb = s(t)r (27) 

Mb = ^e b Rl (28) 
U h = H h R h = s(t) r (29) 

As already mentioned, when the perturbation is well outside the cosmological 
horizon, all of the hydrodynamical quantities have nearly homogeneous profiles with 
their perturbations being small deviations away from the uniform solution. It is useful 
therefore to parametrise the scale of the perturbations using a dimcnsionless parameter 
e that gives explicitly the ratio of the cosmological horizon scale Rr = H^ 1 to the 
physical length scale of the perturbation Rq = s(t)rQ. 



RhV ( 1 V 1 



Ro J \H h sr J s 2 r 2 



« 1 (30) 
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From this equation it is clear that e is a function of time only, and its time derivative, 
as follows from Q24JI. is given by 



-2- + 3(1 +7)- = (I + 37)- 



(31) 



e s e\, s s ' s 

It is important to point out that expression (|31[1 is a consequence of the energy 
conservation law and holds also when 7 is a function of time. 

We treat e as a small parameter for our asymptotic solution and we will see that it 
cancels out in the final equations of this section. However, in numerical computations 
we use an explicit value of e, which for self consistency should be sufficiently small. In 
this case we can keep only first order terms with respect to e in the hydrodynamical 
equations l(T7|) - (|2*2*|l . Then the initial perturbations are defined as: 



R = R b (l + eR) 


(32) 


U = H b R(l + eU) 


(33) 


b- . R ' (l + c6) 


(34) 


a = 1 + ea 


(35) 


e = e b (l + ee) 


(36) 


M = ^ne b R 3 (l + eM) 


(37) 



The tilde quantities are in general functions of both r and t. It is convenient to define 
a new independent time variable £ = ln(s) given by 

is -I- « 38 » 



Starting from equation (fTTjl . and using expressions Q32[l. (|3*3^l and (pT5|> we get 

- + (eRY = -(l + ea)(l + eU) 
s s 

, then linearising this equation and using (|31|l and l|38l) we have 

dR 



(l + 3 7 )i? + 



a + U. 



Similarly, perturbing equation l(TH|) we get 



(eft)" 



a'U 
'IF 



and then, using (|31|l . we have 

^r 9b 
(l + 37)&+^ 



-ra 



Perturbing equation l|19|l and integrating it over r we have 

(1 + 7)a + 7 g = 0. 

Perturbing equation H20JI to the first order in e we get 



-(1 + eM) 



^ + 3§ + (eMY 
e b R 



= - 7 (l + ee)- 



(39) 



(40) 



(41) 



(42) 
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and, using i|24|) . (42) can be re- written as 

- 3(1 + 7)- + 3(1 + 7)- - 3(1 + 7)-eM + (eM)' + 
s R s 

+ 3-eM + 37-ee = 0. 
il it 

Using equation (|2*7jl we get 

eM + eM - 37-eM + 37-eg + 3(1 + j)(eR + eR) = , 
s s 

and finally using equations l|3*T|) . (|3^|l and (|4"T|) we have 

M+^ = -3(l+ 7 )ET. 



(43) 



Perturbing equation l|21(l we obtain 



(l + eM)^ + i e M' = (l + e e)^ 



that can be re-written as 
1 



z = ^M)'. 



Perturbing equation i|22[l . we have 



(1- A'(r)r 2 )(l-2e6)-l = — R 2 (l + 2eU) 



8Tre h R 2 



(44) 



(1 + eM), 



that using (|2*5|) becomes 

- K{r)r 2 + 2eb{l~ K(r)r 2 ) 
and now using expression (|30l) for e we have 



R 2 ^e(2U - M). 



K(r)r 2 + 2eb(l~ K(r)r 2 ) =— (2U - M) . 



Since e«l, one can drop the last term on left-hand side, obtaining finally 
1 



U 



M - K{r)r^ 



(45) 



(46) 



This is the only equation of the system where K(r) appears explicitly. From the 
definition of T given by (|16|l one can rewrite (|46|l as 

1-r 2 

K{r) = —^-. (47) 

This relation shows also the connection between the profile K{r) and T, which is an 
invariantly defined quantity representing an energy per unit mass. 
Substituting l|46|l into equation Q43JI we get 

dM 5 + 3nr,-. 3. „ (4g) 



^Jfr=!(l+7)*(r)r8, 



and from this expression we can see that it is possible to separate the variables (r, £) 
because, as we have assumed at the beginning, 7 is just a function of time. We 
therefore write 



(49) 
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that gives the following differential equation for the function <&(£) 

^ + — *=2<l+7). (50) 

Expression (|49|) is the solution for the mass perturbation M and inserting this into 
14411 , (|41|l and l|46[l we obtain the following expressions for first order perturbations of 
energy density, e, lapse, a and radial velocity, U : 



e = m^[r 3 K(r)]'rl 



U 



MO ~ 1] K{r)rl 



(51) 
(52) 
(53) 



Substituting these equations into and il4TJ|) we find two differential equations for 
R and b, given by 



-Mt)-l]K{r)rl 



^r db 
(l + 3 7 )6+^- 



3^ VK{r))' 



(54) 
(55) 



d£ (1 + 7) 

that can be solved introducing two new functions of time /i(£) and /2(£) and writing 



1 



MO*- 



3r 2 



J X(r))' 



MO 



if(r) 



(56) 
(57) 



The relations between /i,2(£) and <!>(£) are obtained by introducing i|56|) and l|57(l into 
(59 and |nSl>: 



+ (l + 37)Ji(0 



7 



1 + 7 

+ (1 + 37)I 2 (0 = [*(0-1] 



(58) 



(59) 



This completes the system of equations. It worth mentioning that in these equations, 
time evolution is given by e and while all spatial dependencies are given by K (r). 



3. Perturbations in a multifluid medium 

In general the universe consists of different fluid components characterised by different 
constant values 7 j. We can introduce an effective 7 defined as 



7(s) 



(60) 



which, in general, is a function of time. At some moment of time to (which is not 
necessarily the same moment when we impose the initial conditions) we specify the 
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fractional contributions /,; of the different fluid components with the corres 
coefficients 7,, having for each component 

e t o \ s 
= he 

i 

Using expressions flfiTJl and (ffi2|l one can write 

s \ -3(1+7) 

e-i = ho e — 

then from i|ti(J|) and the equation of state © we have 
= E^o7^- 3(1+ ^ 

We now define a new quantity Q 

Q( S )e^/ i0S - 3 ^)«( S ). 

i 

Taking the time derivative of this quantity, one gets 

;f = s f - - 3 D 1 + 5-3(1+70 = - 3 ( x + ^ 

' i 

and finally Q can be related to H b by the following differential equation: 
1 dQ 1 dHfe 

Using expressions (|67l) and (j68(l into H50JI one obtains 
/ _ J_ dflA $ = _ 1 dff 6 



d£ V H b d£_ J H b d£_ 

The solution of this equation can be written in the form 
$ = Fe x , 

where F and [i satisfy the following differential equations 



dF 1 dH, 



which give 



dt H b rfe 

dX 1 dH b 

dt, H b dt, 



-A _ lns-laffi, = 5 
H b J H b ' 



and finally 



s Jo 



Ms) 
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Now we can solve the two differential equations for Ii and I2 that are both of the 
same form, 

s ^s) 2 
ds 

where Fi(s) = jj^y ^(s) and i*2(s) = — 1- From the expression l|31[l we see that 

1 de 00/ d(l/H^s 2 )\ , . 

and this allows the differential equations i|71|) to be written in a form that can easily 
be integrated, 



d 
ds 



' (73) 



therefore 

/( 8 )=^(.),^*F(,) BK L- 5 i, (74) 
which gives the following expressions for the two functions 1\ and ^2 : 

The initial conditions for all physical quantities are now determined in a self 
consistent way: the functions Q(s) and Hb(s) specify the multifluid medium consisting 
of different species of particles and fields (for example massive particles, photons, 
scalar fields, cosmic strings, etc.) which are present in the Universe, and hence in 
the configuration, at any time preceding the time chosen for imposing the initial 
conditions. In the particular case of a single fluid component, with 7 = const, the 
functions Ji, and I2 are also constants. We have in this case 

H(s) cx s-^t 11 , (77) 

and 

5 + 37 

h = 1 $ = ^2 (79) 

(1 + 37XI + 7) (l + 3 7 )(5 + 3 7 ) 1 1 

/ ^TT3^ ((& - 1 ^- (l + 37)(5 + 37) (80) 
Therefore in this simple case the tilde-quantities are time independent and the time 
evolution is determined only by the parameter e. From (|31f) we can calculate the time 
evolution of e 

2(1 + 3 T ) 

e(t) cx ( -J , (81) 

which is the same as the standard solution for a growing mode in cosmological 
perturbation theory (see for example |31M34| ). 
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4. General properties of the curvature profile 



4-1- The curvature profile and the comoving curvature perturbation 

The curvature profile K(r) used in the present paper is linked to the comoving 
curvature perturbation TZ frequently used in literature (see for example, HJ). The 
simplest definition of TZ comes from the theory of linear perturbations in the FRW 
metric, where the spatial part of the metric tensor is given by 

9ij = s 2 [(1 + 211) Sij] . (82) 

The calculation of the spatial scalar curvature gives 



i?(3) = 4 *!ft 



(83) 



where k is the comoving wavenumber associated with the perturbation. It is possible 
to generalise the definition of TZ for non linear perturbations of the metric |32l l3~3"| , 
with TZ being still time independent (as first shown in 1963 by Lifshitz & Khalatnikov 
|30j ) and converging to H82|) for small perturbations. Starting from the expression for 
the scalar curvature given by (|83l) for generic metric perturbations one can see that 
the comoving curvature perturbation TZ is connected with the pure growing mode in 
the case of small density perturbations 

Se 2(1 + 7) / k \ 2 

eb 5 + 37 \sH J 
where 7 is constant. This relation is obtained using a quantity called the "peculiar 
gravitational potential" <I>g defined as 

Se _ 2 / k x :> 
eb 3 \sH y 
and linked to TZ by a first order differential equation 



(85) 



1 



5 + 3 7 



-(1+7)^. 



(86) 



H 2 

This equation has a pure growing solution given by 

*o = -¥r?^ (^) 
5 + 37 

where TZ and $0 are time independent when the length-scale of a perturbation is 
much larger than the horizon scale. If we calculate the spatial curvature scalar (see 
appendix) in terms of the curvature profile K(r), we find that 



R = 6 



37-I 
: 5 + 3 7 



K. 



2(2 + 3 7 ) ^ 



K 

„2 



where 



K = K(r) + -K'(r) . 



For the diagonal metric Q the three curvature is given by 

6 r, 



R& 



K(r) + -K'(r) 



(89) 



(90) 



because only the last term in the square brackets of equation (|88|l corresponds to the 
diagonal spatial component of the Ricci tensor. We now identify 

fc = l/r , (91) 
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which is obvious to do because the wave number represents the inverse length of the 
perturbation. Comparing equations 1(83(1 and (|90|l . one can then see that the two 
equations are the same if 

r 



17 ^2 

Tt — —T< 



K(r) + -K'(r)\ . (92) 



2° 

From !jUT)l and (|3T)|l it follows also that 

(93) 



k 

7h / 

and relation 1(92(1 is obtained again by comparison of ((51(1 with 1)84(1 . We can notice 
also that equation ((86(1 is the same as (15011 since the peculiar gravitational potential 
defined in ((85(1 can be written in our terminology as $g = — (2/3)e and we can then 
use equation ((44(1 to express the variable M in terms of e in l(5U(l . We have therefore 
shown that 1Z is directly linked to K(r). 

4-2. Physical properties of the curvature profile 

From the definition of b given by 1(34(1 we obtain obvious mathematical requirement 



1 



1 - K[r)r A > => K{r) < . (94) 

This corresponds to the physical condition that a perturbed spherical region of 
comoving radius r should not be causally disconnected from the rest of the Universe. 
Another important requirement is related to conservation of total mass, given by 



4irr 2 e(r)dr = , (95) 



that implies 



lim r 3 K(r) =0. (96) 



This means that the Universe remains spatially flat at infinity. As explained in 6 , if 
condition ([?]) is not satisfied one unavoidably confronts a violation of the causality 
principle. 

An important parameter characterising the curvature profile is the value of ro 
that specifies the comoving length-scale of the overdense region in the configuration. 
From 1(51(1 we have 

e(r ) = => K(r Q ) + jK'(r ) = 0. (97) 

It is useful to introduce an integrated quantity 6 that measures the mass-excess inside 
the overdense region, as frequently done in literature. From equation j£lj, one can 
see that any spherical integral of e does not depend on the particular curvature profile 
used, but depends only on the value of K(r) at the outer edge of the configuration. 
The expression for 5 is then given by 



5 = -tt4 / 4tt -r z dr = e(s)<P(s)K(r )r{; . (98) 

V 3 / Jo e b 

We can see from this expression that 5 is a gauge invariant quantity since it is directly 

proportional to if(ro)r§, which is also a gauge invariant quantity. Note once more 

the perfect separation of the different physical aspects of the problem: e specifies the 

evolution of the perturbation amplitude in the long wavelength limit, $ represents the 

equation of state and K(r)r 2 specifies the spatial profile of the configuration. 
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5. Parametrisation of the curvature profile 

In this section we specify the profile K (r) in terms of two free parameters. First 
we want to point out that both r and K are gauge-dependent quantities. Therefore 
we need to be careful when transferring the initial perturbation profiles onto the grid 
used for the numerical calculations. The circumference coordinate Rb used in the code 
is gauge invariant and is related to r by Rb = s(t)r, while particular values of s(t) 
and r are gauge dependent. Another gauge invariant quantity is the product K(r)r 2 . 
These two gauge invariant products are connected by the presence of the comoving 
coordinate in both of them and therefore a particular specification of one of these 
three variables will specify the others. 

To do this we compare the length scale Rq of the configuration with the 
cosmological horizon Rh 

Ro = NR H , (99) 

where N is the number of horizon scales inside the configuration. From the definition 
of e in (|30|l , we have 

e =^2' ( 10 °) 

To have e <C 1 we need to have N 3> 1 and N will be one of the input parameters 
used to impose the initial conditions. Expressions i|99|) and (|30() relate the scale factor 
with r Q 

•(t) = ^, (101) 
and we can use this to define the comoving coordinate as 

r = jRb = _^o_ ( } 

s{t) NR H V ' 

where the value of Rb is calculated from equation (J3J 

* b =(i^r< d03) 

and ro is determined by condition (|97|l once the curvature profile is specified. 

To parametrise K(r) we should decide what is a reasonable perturbation shape. 
It is natural to start with a centrally peaked profile of the energy density, such 
that outside the overdense region there is an under-dense region which becomes 
asymptotically flat. To obtain an energy density profile with these properties, one 
should choose a centrally peaked curvature profile which is described in terms of a 
suitable continuous function and tends asymptotically to zero as r — > oo. Continuity 
should be ensured at least in the first and second derivatives, because these derivatives 
play a crucial role in the expressions for the perturbation profiles of density and 
velocity. To preserve the standard normalisation used in cosmology where a closed 
universe has a curvature equal to 1, we choose to normalise the curvature profile to 
have K(0) = 1. 

We start with a family of curvature profiles based on a Gaussian shape, described 

as 

„2 



K(r)= [l+a— lexp -— ) (104) 
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where a and A are two independent parameters. In the particular case a 
profile K(r) is exactly Gaussian. The first derivative of K(r) is given by 



K'{r) 



r 

A 2 



a- 1 



"A 2 



exp 



2A 2 



0, the 



(105) 



and we use this to calculate the value of r$ as function of a and A using l|104|l and 
H97fl . We obtain a quadratic equation with the following solution for ro: 

rl = 3A 2 if a = (106) 



2 (5a 



2) + y/{5a - 2) 2 + 24a 



2 a 



if a ^ 



(107) 



that we substitute into (|102f> to express the spatial comoving coordinate r in terms 
of the background coordinate i?b- It is then convenient to express ro in terms of the 
initial perturbation length scale i?o- We therefore introduce a dimensionless ratio of 
the circumferential radius to Ro defined as Z = Rb/Ro, which is equal to 1 at the edge 
of the overdense region. Thus the exponential argument in (102) is given by 

r 2 _ F(a) 
A 2 ~ 

where 



2a 



-Z 2 



(108) 



F(a) = (5a - 2) + V '(5a - 2) 2 + 24a if a ^ 



(109) 



F(a) 



3 if 







(110) 



and one can see that the parameter A cancels in the expression for r/ro, making 
clear the physical meaning of a and A: for the energy density perturbation profile 
expressed in terms of Z we get 

2 



1 



5 



i 



2a 



Z 2 - r: 



a f F(a) 
2 I 



2a 



x exp 



F(a) 
4a 



Z' 



(111) 



and we see that the spatial profile given by the expression inside the parentheses 
depends only on a. Hence A, appearing only outside the brackets, parameterises the 
amplitude of the density perturbation. The separation of factors in (|111|) allows us to 
control the values of the input parameters independently. 

The value of the second derivative of K(r) at the centre is given by 

K "(°) = ^T- ( 112 ) 

When a > 1, the value K"(0) is positive and an off-centred peak appears in the 
profiles of both curvature and energy density. For a < 0, there is a second overdense 
region, corresponding to a second solution for r - Because we are interested in solutions 
that represent perturbations centred at r = surrounded by an underdense region 
which becomes asymptotically flat, we restrict our attention to the range of shapes 
given by < a < 1. 

In Figure ^ we plot corresponding curvature and energy density profiles for three 
different values of a between and 1, keeping the same value of A. 
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2.5 



Figure 1. The left hand plot shows the curvature profile K (r) as a function 
of the comoving coordinate for three different value of a (0, 0.5 and 1). 
The right hand plot shows the corresponding profiles of the energy density 
perturbation e plotted as functions of Z. These cases and those described 
in the next figures have been calculated with 7 = 1/3 (equivalent to 
$ = 2/3.) 



The particular case with a = in (|lllfl coincides with one of the profiles used in 



e(Z) 



2 



,xp(-|z 2 ). 1118) 



where the parameter A has the same meaning as the quantity used in |19l I25| to 
parametrise the amplitude of the density perturbation. 

An advantage of the present description is the possibility to calculate the 
perturbation amplitude S, defined in (|98|) . directly from the curvature profile. Using 
l|113|) for an arbitrary value of a we have 

1 A 2 , J F(a)\ ( F(a)\ , , 

Another important point is given by condition 194JI , which should be taken into account 
to calculate the maximum allowed value of the density perturbation amplitude for each 
profile. Therefore one needs to calculate the maximum value of A for each a. The 
results of these calculations are presented in Figure [2] for two different cases (a = 
and a — 1). From these plots, we can see clearly how the density perturbation 
amplitude S varies with A. 

Next we consider another parametrisation of curvature that allows us to work 
with curvature profiles with a sharper transition from K = 1 to K = 0. The new 
parametrisation is given by the following expression: 

1 if r<A» 



1+ 2A2 j eXP 2A^I if r>A * 
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Figure 2. The left hand plots show the curvature profiles K(r) as 
functions of the comoving coordinate r, while the right hand plots show 
the corresponding profiles for the energy density perturbations e plotted as 
functions of Z. For a = (upper plots) the different profiles correspond 
to A between 0.2 and 1.16, which is the maximum value allowed by (19411 . 
The values of A are higher for the higher curves. For a = 1 (lower plots), 
the values of A used are between 0.15 and 0.77. 



and its first derivative is 

if r < A* 

K ' {r) = < (r - A») 3 / ( r -A.n . f . ^ 

In figure |21 we present the parametrisation given by (|115J) . indicating explicitly the 
two independent parameters A and A*: A* specifies the radius of the plateau with 
K(r) = 1, A describes the sharpness of the transition region. In particular, when 
A* = the curvature profile coincides with the previous case for a = 1. Applying 
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Figure 3. This figure shows the behaviour of K(r) given by 1115H . 
indicating expiicitly the meaning of the two parameters A and A*. 



condition Q97JI we obtain a quadratic equation to calculate the value of ro, that we 
solve numerically using the Newton- Raphson method |35| . The maximum value of A* 
satisfying condition l|94|l is 1, while the minimum value of S, obtained when A = 0, is 



In Figure 0J the top left plot shows the dependence of K(r) on A with A* kept 
fixed, while the bottom left plot shows the dependence of K(r) on A* with A kept 
fixed. In the right hand plots, we show the corresponding profiles of e(r). 

6. Numerical tests and calculations 

6.1. Numerical tests 

Next we discuss the numerical tests performed so as to be confident that the 
perturbation profile K(r) has been introduced consistently into the code. For this we 
analysed the accuracy of our first order approach (first order with respect to e — /N 2 ). 
In other words, imposing a perturbation with a length scale much larger then the 
cosmological horizon r# , we checked whether e is small enough to make higher order 
terms negligible. A first requirement for the necessary value of N is that the comoving 
lengthscale of the horizon r# should be sufficiently smaller than the sharpness A of 
the curvature profile. This means that rn << A. From equation (|102fl . r# = ro/N 
which then gives the condition, 



Using the curvature profile given by l|l(J4|) we have ro — A for < a < 1, while 
using the curvature profile given by l|115|) we have that ro ~ A* + A. Substituting 



5, 



2A 2 J3N 2 . 




(117) 
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Figure 4. The left hand plots show the curvature profile K(r) for A* 
constant and varying A at the top row and for A constant and varying 
A* at the bottom. The corresponding behaviour of the energy density 
perturbation e is shown in the right hand plots. 



these two relations into (|117fl . we see that both conditions are satisfied if 

N » 1 + . (118) 

This expression is valid also for the curvature profile given by i| 104(1 if we put A» = 0. 
From this expression we can see that a sharper profile, with A„/A 3> 1, requires 
larger values of N. 

For each set of initial input parameters we apply several numerical tests. Here 
we present four sample cases that demonstrate how these tests work. The main test 
is based on equation (|47() . which is used to rederive K(r) after calculating the initial 
conditions for the various quantities. Expression (|47() is valid when the higher order 
terms are negligible and so we expect to find some difference between the profiles of 



Curvature profiles as initial conditions for primordial black hole formation 



20 



1 Efi 



0.75 - VJ 



a= 



0.25 - 



0.75 



a= 1 



0.25 




Figure 5. These plots show the comparison between the analytical (solid 
line) and the numerical (dashed line) profiles of curvature given by I104H . 
The different dashed lines correspond to TV = 1,2,3,4,6,8,10 with the 
higher curves corresponding to the higher values of TV. The profiles in the 
left hand plot are characterised by (a = 0, A = 1.15), and the profiles in 
the right hand plot by (a = 1, A = 0.77). 




Figure 6. These plots show the comparison between the analytical (solid 
line) and the numerical (dashed lines) profiles of curvature given by 111511 . 
The different dashed lines correspond to TV = 1, 2, 4, 5, 7, 10 in the left plot 
and TV = 1, 5, 25 in the right one, with the higher curves corresponding to 
the higher values of TV. The profiles in the left hand plot are characterised 
by A*/A = 1, and the profiles in the right hand plot by A*/A = 20. 



K(r) calculated with (|47|l and the original ones introduced analytically to describe 
the initial conditions. In Figure [S] we plot the result of this test using K(r) given by 
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H X U4|) . with a — on the left and a = 1 on the right. In FigureEJwe plot K(r) given 
by (|115fl . using the ratio A* /A equal to 1 in the left panel, and equal to 20 in the 
right panel. In all of these plots the solid line shows the analytical profile of K(r), 
while the dashed lines corresponds to the profiles of K(r) calculated with l|47|) using 
different values of N. The first value used was N = 1, that we know is not large 
enough, and this is also clear from the plots; we then increased N. In the first three 
cases, corresponding to curvature profiles with A* /A < 1, we found that N = 10 
is large enough to achieve acceptable precision. In the last case, when A* /A = 20, 
we found that N = 25 is sufficiently large. Because we did not see any significant 
difference between the analytical and numerical values of K(r), for the highest values 
of N plotted in Figure [5] and we assume that these values of N are satisfactory. To 
be more certain of this we also made another test, analysing numerically the evolution 
of configurations corresponding to the curvature profiles K (r) presented in Figure 03 
and EJ We have verified that the values of N indicated by the previous test do give 
essentially the same numerical evolution as for larger values of N (a standard test for 
numerical computations). 

To complete our tests we used the constraint equation l|16l) to estimate the 
precision of the code and we found that, with the number of gridpoints that we used 
in our simulations, the numerical scheme has accuracy better than one part in 10 4 (by 
the numerical scheme here we mean the code itself plus the initial conditions.) 

6.2. Description of the calculations 

The calculations of black hole formation performed with initial conditions specified 
with a curvature profile, show similar hydrodynamical features as those presented in 
25 , where the same numerical technique was used. We present here four sample 
cases: two correspond to black hole formation, where the observer time slicing has 
been used, and two correspond to no black hole formation, where we have used the 
cosmic time slicing. Each case is characterised by a particular value of the density 
perturbation amplitude 8 and it is useful to introduce at this stage a measure of 5 
which is independent of the initial scale used when imposing the initial conditions. 
Looking at expression 1|98[) . we can see that in the case of constant 7, if we define 
5 = N 2 5, the expression becomes time independent, giving in the radiation epoch 

5 = |jf(r )rg ■ (119) 

Comparing this value with the value of 6 at horizon crossing, we find that it is of 
the same order (with a small correction due to the non linear growth of the energy 
density near to horizon crossing) . Therefore for the rest of the discussion we use S as 
the amplitude of density perturbations. The dividing line between perturbations that 
collapse to form black holes and perturbations that disperse into the background is 
characterised by a threshold amplitude 6 C . This value was estimated numerically by 
making a converging sequence of numerical calculations until a satisfactory precision 
in determination of 8 C was achieved. 

The observer time coordinate u is defined by 

/ du = adt - bdr (120) 

where / is the lapse in the new coordinate system. The metric is then 

ds 2 = -f 2 du 2 - 2fbdrdu + R 2 {d0 2 + sin 2 8d(j) 2 ). (121) 
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Figure 7. A typical evolution leading to black hole formation: the initial 
curvature profile used is characterised by a = and A = 1.02, which 
give (S - So) = 1.3 x fO" 2 . The upper plots show the profile of radial 
velocity U/c and lapse / at different times, while the bottom panels show 
the corresponding profiles of 2M/R and mass M. 



When using observer-time slicing, black hole formation occurs when both the lapse / 
and r + U tend to zero (see US]). These conditions are reached asymptotically and 
correspond to the formation of the event horizon, seen by a distant observer with an 
infinite redshift. 

We illustrate the two examples of black hole formation, shown in Figures [7| and 
00 by plotting profiles of the radial velocity U, lapse /, 2M/R ratio and mass M as 
functions of the circumferential radius R at different times. These show the formation 
of black holes for different values of (6 — 6 C ) with the smaller value giving rise to the 
lower mass black hole as expected [H| ^3 EH Ei] . 
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Figure 8. A typical evolution leading to black hole formation: the initial 
curvature profile used is characterised by A* = 0.8 and A = 0.f8, which 
give (5 — 8 C ) = f .8 x fCF 3 . The upper plots show the profile of radial 
velocity U/c and lapse / at different times, while the bottom panels show 
the corresponding profiles of 2M/R and mass M. 



The upper left hand plots of these two figures show the behaviour of the radial 
velocity U. Here the time sequence starts from the top curve and continues towards 
the bottom one. Note that, after a certain time appreciate how U becomes negative 
in the central region, corresponding to the region where the black hole forms. In 
particular the event horizon, as can be seen clearly looking at the plots for 2M/R, 
forms at the minimum of the radial velocity profile. Outside this region there is first 
an intermediate shell with negative velocity corresponding to matter accreting into 
the black hole, then there is a region of positive velocity corresponding to the rest of 
the expanding universe. 
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The upper right hand plots show the behaviour of the lapse / with the time 
sequence of the curves given by the decreasing central value of /. One can sees that 
a plateau forms in the central region, where the black hole is forming, and that the 
lapse tends asymptotically to zero there. The "freezing" of the evolution, as mentioned 
above, is one of the fundamental criteria for black hole formation. The other important 
feature to ensure us about black hole formation is the trapped surface condition 
(r + U -> 0) which corresponds to 2M/R -> 1. In the left bottom plots 2M/R is 
shown asymptotically tending to 1 at the maximum of the curve. This determines 
the location where the event horizon will form. In this case the time sequence of the 
curves is given at first by the decreasing values of 2M/R at the right edge of the plot, 
and then by the increasing of the maximum of these curves towards 1. 

Finally the bottom right hand plots show the corresponding behaviour of the 
mass M, indicating the value of the mass of the black hole. At late times, one sees the 
convergence of the curves. In figure [7| and [H] the circumferential radius R and mass M 
have been normalized with respect to the values for the cosmological horizon scale at 
the horizon crossing time. 

Results from two representative calculation with S < 5 C (using cosmic time) 
are shown in Figure where we have plotted the energy density normalised with 
respect to the current background value, as a function of R/Rh(U). The left hand 
plot corresponds to a small curvature perturbation far from the threshold 6 C with 
S = 4.46 • 1CP 3 and (6 — 5 C ) — —0.45. In this case negative velocities are never 
reached: initially the slower expansion of the central region proceeds more slowly than 
the expansion of the background. After horizon crossing, the perturbation is damped 
and compression a wave propagates outwards through the expanding medium. The 
other example shown in the right hand plot shows a case with [6 — S c ) = —2.67 • 10 -2 . 
In this case the perturbation amplitude 5 = 0.49 is close enough to S c = 0.52 
to give rise to an initial collapse of the central region. Because the gravitational 
potential is not large enough to form a trapped surface, as happens when 6 > S c , the 
collapsing matter bounces producing a compression wave which expands outwards. 
These hydrodynamical features were first seen in |S] and were then analysed in detail 
in 

As a final result of these computations we determined the parameter ranges 
corresponding to black hole formation for the curvature profiles given by (|1C)4|) and 
H115J) . Figure El summarises the results. The left hand plot shows the results obtained 
for the value of A and a used in parametrisation (1 104(1 , with the plane being divided 
into three different regions: the first one corresponds to no black hole formation, the 
second one to black hole formation and the third one to disconnected configurations 
(K(r)r 2 > 1). It can be seen that the region of parameters corresponding to black-hole 
formation is quite narrow. 

The same results are presented in terms of S and a in the right hand plot of 
Figure ITU1 This shows very weak dependence of S c and 5 max on a. We found that 
S c ~ 0.45, <5 max ~ 0.60 for a = and S c ~ 0.47, S max ~ 0.62 for a = 1. 

In Figure ITTI we present similar results for the second curvature parametrisation 
given by 1(115(1 . The left hand plot shows the parameter space for A and A*, and one 
can see the same three regions in the left hand plot as in Figure El The right hand 
plot shows the same results in terms of A* and 8. 

For wide range of parameters explored so far, we have not yet seen formation 
of shocks as had seen in some earlier works jHJ El CHO EH ■ Work is in progress to 
investigate this further. 
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Figure 9. The left hand plot shows a typical evolution of the energy density 
profile of a perturbation with 5 <C S c . In this case we have used expression 
itToH with a = 0, A = 0.1 and TV = 10, corresponding to 8 = 4.6 • 10~ 3 . 
The right hand plot show the typical evolution of a collapse when S < S c . 
In this case we have used expression 11151 with A, = 0.5, A = 0.37 and 
N = 10, corresponding to S = 0.49. 
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Figure 10. These plots show which values of a, A and <5 lead to black hole 
formation or to an initial perturbation already disconnected from the rest 
of the Universe for the case of K(r) given by expression 110411 . 



7. Conclusion 

Using the quasi homogeneous solution, we have impose initial conditions for PBH 
formation by introducing a time independent curvature profile K(r), related to the 
curvature perturbation 1Z used elsewhere in the literature. This profile is the only 
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source of perturbations in the fluid quantities and those perturbations behave as pure 
growing modes. We have obtained an analytical solution of the Misner - Sharp system 
of equations expressing perturbations of density and velocity in terms of the curvature 
profile K(r) in the rather general case where matter in the Universe can be treated as 
an arbitrary mixture of perfect fluids. 

We have performed numerical calculations for 7 = 1/3 to test the self consistency 
of these initial conditions for two different parameterisations of K{r). 

Specifying the initial conditions using a curvature profile seems to be an 
improvement in the analysis of the PBH formation scenario since it allows initial 
conditions for all the hydrodynamical variables relevant for the problem to be specified 
in self consistent way. 

We have shown that, in agreement with [2|, the formation of PBHs requires higher 
amplitudes of strong metric deviation for sharper profiles of K(r). Thus the Gaussian 
curvature profile, which is the least sharp, is the most favourable for PBH formation. 
We found that in this case S c ~ 0.45 in agreement with [21]. The comparison of 
results obtained for two different parametrisation of curvature profiles has shown 
that the shape of the initial profile plays crucial role for PBH formation. This is 
an important reason to calculate the probability for having different curvature profiles 
(work in progress) [37] and to link this probability with statistics of PBHs. Such a 
link should greatly improve the constraints on possible cosmological models obtained 
from observational upper limits on PBHs. 
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Appendix A. Scalar curvature 

The general spherically symmetric metric in the cosmic time coordinate is 

ds 2 = -a 2 dt 2 + b 2 dr 2 + R 2 dn 2 (A.l) 

where the metric tensor g a p is diagonal and a, b and R are, in general, functions of r 
and t. The Christoffel symbols are calculated by 

( dgp a dg l5 dg 01 \ 

while the components of the Ricci Tensor are given by 

<5r 7 r)r 7 

it a p - L Sy L a p-L S0 1 aj + ^-g- (A.6) 

and to calculate the Ricci scalar we need to take into account only the diagonal 
components R aa . 

In the background universe with the FRW metric the metric components are 
therefore 

s(t) 

a = l b=-j=LL= R = s(t)r, (A.4) 

where the curvature if is a constant equal to 0, ±1. The diagonal components of the 
Ricci tensor are 

i?oo = 3- (A.5) 
s 

_ ss + 2s 2 + 2K 
Rl1 ~ l-Kr* (A - 6) 

R 22 = (as + 2s + 2K ) r 2 (A.7) 

R33 = #22 sin 2 *? (A.8) 
and these are used to calculate the scalar curvature 

*-•(; + ? + ?)■ <"> 

While the three curvature is obtained by removing the terms with the time derivatives, 

R (3) = e4- (A.10) 

s 1 

In the general case of the metric given by (| A. 1|) . the diagonal components of the 
Ricci tensor are 

1 da / 1 db 2 dR\ a da ( 2 OR 1 db 
00 ~ a di (bdi + R~dt ) + tf'dr' \R~dr~ ~ bdr 



a d 2 a 1 d 2 b 2 d 2 R 

+ ]pq^2 - lap - -j^-Q-p 

b db ( 2 OR lda\ 1 db (I da 2 dR 
' 1 a 2 dt \R dt a dt j b dr \a dr R dr 



(All) 
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R-2 



RdR fldb 
a 2 ~dt [bdi 

ldR 



R: 



a dt 
R22 sin 2 



1 da 
a dt 

ldR 
b~dr~ 



b d 2 b 1 d 2 a 2 d 2 R 
a 2 dt 2 a dr 2 R dr 2 

RdR fldb Ida 
b 2 dr V b dr a dr 



(A.12) 



R 



1 d 2 R 1 dR 

a 2 dt 2 b 2 dr 2 



To write the scalar curvature it is useful to define the following operators: 



D, 



ld_ 

a dt 



D 



1_ (P_ 



D r 



18_ 

b dr 



Dt 



b 2 dr 2 



1 (A.13) 
(A.14) 

(A.15) 



We then get 



R 



D 2 b D 2 R 
— + 2^— 
b R 



D 2 a 2 D 2 R 



2 £taDtR _ DtaDtb 
a R a b 



R 

D r b D r R 
~b ^~ 
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D r a D r b 



D r R 
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DtbDtR 
' b R 



D r a D r R 
'~~a ^~ 



1 

R 2 



(A.16) 



and substituting the expressions for the coefficients a, b and R given by H32I34I35J1 and 
expressing a,b,R in terms of K(r), we get with the use also of (IHlll . that the Ricci 
scalar is given by 



R 



37- 1 
: 5 + 3 7 ' 



K 



2(2 + 3 7 ) ^ 
5 + 3 7 



K_ 

72 



where 



K = K(r) + -K'(r) . 



(A.17) 
(A.18) 



It is useful to compare (|A.17f) with (|A.9|1 . 
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